Analysis of ligand-receptor interactions was done with the CellPhoneDB package, Efremova et al, 2019. The basic analysis was done in Python as described in the paper, whereafter the data were edited and plotted in R. In case of too many ligand_receptor pairs selection of relevant pairs was made. As in CellPhoneDB some pairs are given as receptor-ligand, this was manually transformed to ligand_receptor by computation in excel. The order of clusters was edited to be the same as in other figures.
Python and CellphoneDB were installed as system resources, please see code.
From the original seurat file a raw count file is generated with ensemble symbol annotations and a metadata file with the cluster annotations. Start is: seu_2021_May.R generated and saved in CMC2.Rmd
#install.packages(c("readxl","writexl"))
suppressPackageStartupMessages({
library(igraph)
library(ggplot2)
library(biomaRt)
library(Seurat)
library(dplyr)
library(data.table)
library(reshape2)
library(openxlsx)
library(readxl)
library(xlsx)
library(reshape2)
library(data.table)
#library(reticulate)
library(writexl)
})
path1 <-'/Users/ChristerSylven/Human-Prenatal-Cardiomyocytes/'
path2 <-'/Users/ChristerSylven/'
seu <- readRDS(paste0(path1, 'seu_2021_May.R'))
#seu
#names(seu[[]])
#table(seu@meta.data$sub_cluster_SAN2) # SAN should be n=8Idents(object = seu) <- "sub_cluster_SAN2"
seu_Phone <- subset(x = seu, idents = c('M0','M1','M2','M3','M4', 'M5', 'M6', 'M7', 'SAN','0','2','3','4','6','7','9','10','12','14','15','16'))
seu_Phone # all cell clusters except erythocytes## An object of class Seurat
## 30646 features across 3371 samples within 2 assays
## Active assay: SCT (15323 features, 3000 variable features)
## 1 other assay present: RNA
## 2 dimensional reductions calculated: pca, umap
#names(seu_Phone@meta.data)
seu_Phone@meta.data$sub_cluster_SAN2 <- factor(seu_Phone@meta.data$sub_cluster_SAN2)
table(seu_Phone@meta.data$sub_cluster_SAN2)##
## M0 M1 M2 M3 M4 M5 M6 M7 SAN 0 2 3 4 6 7 9 10 12 14 15
## 190 141 103 79 74 51 57 23 8 652 475 369 340 160 148 122 119 111 76 44
## 16
## 29
#count_raw <- seu@raw.data[,data_object@cell.names]
count_raw <- as.data.frame(GetAssayData(object = seu_Phone, slot = 'counts'))
count_norm <- apply(count_raw, 2, function(x) (x/sum(x))*10000)
counts <- count_norm
ensembl = useMart("ensembl",dataset = "hsapiens_gene_ensembl")
genes <- getBM(filters = 'hgnc_symbol',
attributes = c('ensembl_gene_id','hgnc_symbol'),
values = rownames(counts),
mart = ensembl)
counts <- counts[rownames(counts) %in% genes$hgnc_symbol,]
counts <- tibble::rownames_to_column(
as.data.frame(counts), var = 'hgnc_symbol')
counts <- plyr::join(counts, genes)
counts$hgnc_symbol <- NULL
counts <- cbind(counts[,which(colnames(counts) == 'ensembl_gene_id')], counts)
colnames(counts)[1] <- 'Gene'
counts$ensembl_gene_id <- NULL
metadata <- data.frame(Cell = rownames(seu_Phone@meta.data),
cell_type = Idents(object = seu_Phone))
write.table(counts,
file = './out/counts.txt',
quote = F,
col.names = T,
row.names = F,
sep = '\t')
write.table(metadata,
file = './out/metadata.txt',
quote = F,
col.names = T,
row.names = F,
sep = '\t')Normalized counts (ensembl_gene_id) and metadata tables are formatted according to Cellphone DB CellphoneDB is then run and heatmap for all ligand_receptor interactions between all clusters generated in the user/out folder:
#library(reticulate)
Sys.setenv(RETICULATE_PYTHON = "/Users/ChristerSylven/Library/Python/3.8/bin/")
# install cellphonedb as a system resource and use path as below
system('python3 -m venv cpdb')
system('source cpdb/bin/activate')
system('pip install cellphonedb')
# to generate means and p values and heatmap between all celltypes use:
# system("cellphonedb method statistical_analysis metadata.txt counts.txt --counts-data=gene_name --iterations=100 --threads=2 --result-precision 2 --output-path=/Users/ChristerSylven/out/")
# system("cellphonedb plot heatmap_plot metadata.txt --pvalues-path=./out/pvalues.txt")
# alternative: use cellphone call without path directly in terminal after having installed cellphonedbTo show all ligand_receptor interactions between of non-CM cell types of relevance according to scRNAseq deconvolution on ST maps in OFT (M7) and atrial Conduit (M5) tt containing spots, metafiles for these cell types are saved
#heatmap plot and count_network (table for heatmap) for subgroups
# Metadata for groups
metadata <- read.table('metadata.txt')
colnames(metadata) <- metadata[1,]
metadata <- metadata[-1,]
cc <- c('M7', '2', '7', '6', '0', '3', '15') # OFT v2
metadata_M7 <- metadata[metadata$cell_type %in% cc,]
metadata_M7[,2] <- factor(metadata_M7[,2]) # delete cluster w zero cells
write.table(metadata_M7,
file = './out/metadata_M7.txt',
quote = F,
# col.names = T,
row.names = F,
sep = '\t')
metadata <- read.table('./out/metadata_M7.txt')
table(metadata_M7[,2])##
## 0 15 2 3 6 7 M7
## 652 44 475 369 160 148 23
system("cellphonedb plot heatmap_plot ./out/metadata_M7.txt --pvalues-path=./out/pvalues.txt")Following the system command
the OFT heatmap and its count_network data are saved in the out folder
to construct a triangular heatmap of interactions:
# heatmap with desired order of clusters
count_network <- read.table('./out/count_network.txt', header = TRUE)
#OFT cc <- c('M7', '2', '7', '6', '0', '3', '15')
count_wide <- setDT(count_network %>% arrange(factor(SOURCE, levels = c('M7', '2', '7', '6', '0', '3', '15'))))
count_wide <- dcast(count_wide, factor(SOURCE, levels = unique(c('15','3', '0', '6', '7', '2', 'M7'))) ~ factor(TARGET, levels = unique(c('15','3', '0', '6', '7', '2', 'M7'))), value.var = 'count')
count_wide <- count_wide[,-1]
count_wide[lower.tri(-count_wide)] = NA
cc <- c('OFT', 'Epicardium_d', 'Fb_smooth_m', 'Fb_l_vascular', 'Endothelium', 'Fb_skeleton', 'Schwann')
colnames(count_wide) <- rev(cc)
count_wide$celltype <- rev(cc)
count_wide2 <- count_wide
triangular_plot <- function (count_wide2)
{
count_wide2 <- na.omit(melt(count_wide2, 'celltype', variable='cohort'))
count_wide2$cohort <- factor(count_wide2$cohort)
count_wide2$celltype <- factor(count_wide2$celltype)
#head(count_wide2)
count_wide2$log <- log(count_wide2$value+1)
ggplot(data = count_wide2) +
aes(x = cohort, y = celltype) +
ggtitle('') +
theme_bw() +
xlab('') +
ylab('') +
geom_tile(aes(fill = log), colour = "transparent") +
scale_fill_gradientn(colours = c("dodgerblue4", 'peachpuff', 'deeppink4')) +
scale_x_discrete(limits = cc,
labels = cc) +
scale_y_discrete(limits = cc,
labels = cc) +
theme(
axis.text.x = element_blank(),
axis.text.y=element_text(size=12),
axis.ticks=element_blank(),
axis.line=element_blank(),
axis.title.x=element_blank(),
axis.title.y=element_blank(),
panel.border=element_blank()) +
annotate("text", x = (1:length(levels(count_wide2$cohort))), y = 1:length(levels(count_wide2$cohort))- 0.75, label = rev(levels(count_wide2$cohort)), size=4, hjust = 1, angle = 90) +
coord_cartesian(ylim = c(8, -2), clip = "off")
}
triangular_plot(count_wide)In the same way, create a metadata file for atrial conduit and colocalized cell types and generate count_network and triangular heatmap.
metadata <- read.table('metadata.txt')
colnames(metadata) <- metadata[1,]
metadata <- metadata[-1,]
cc <- c('M5', '0', '2', '6', '7', '12', '15')
metadata_M5<- metadata[metadata$cell_type %in% cc,]
metadata_M5[,2] <- factor(metadata_M5 [,2]) # delete cluster w zero cells
write.table(metadata_M5,
file = './out/metadata_M5.txt',
quote = F,
#col.names = T,
row.names = F,
sep = '\t')
table(metadata_M5[,2])##
## 0 12 15 2 6 7 M5
## 652 111 44 475 160 148 51
system("cellphonedb plot heatmap_plot ./out/metadata_M5.txt --pvalues-path=./out/pvalues.txt")count_network <- read.table('./out/count_network.txt', header = TRUE)
count_network[,1] <- as.factor(count_network[,1] )
count_network[,2] <- as.factor(count_network[,2])
count_network$log_count <- log(count_network$count +1 )
count_network$log_count2 <- count_network$log_count
#count_network
#cc <- c('M5', '6', '2', '3', '12', '4', '7')
# to construct a triangular heatmap of interactions do:
#Conduit + SAN
count_wide <- setDT(count_network %>% arrange(factor(SOURCE, levels = cc)))
count_wide <- dcast(count_wide, factor(SOURCE, levels = unique(c('15','12', '7', '6', '2', '0', 'M5'))) ~ factor(TARGET, levels = unique(c('15','12', '7', '6', '2', '0', 'M5'))), value.var = 'count')
count_wide <- count_wide[,-1]
count_wide[lower.tri(-count_wide)] = NA
#count_wide
cc <- c('Conduit', 'Endothelium', 'Epicardium_d', 'Fb_l_vascular','Fb_smooth_m', 'Epicardium', 'Schwann')
colnames(count_wide) <- rev(cc)
count_wide$celltype <- rev(cc)
count_wide2 <- count_wide
triangular_plot(count_wide)The basic /out/means.txt & /out/pvalues.text are first pruned by deletion of duplicates
Then the relation between integrin and non-integrin ligand_receptor pairs is plotted for interaction between different cardiomyocytes and cardiomyocytes and non-cardiomyocytes.
#subgroup output
means <- read.delim("/Users/ChristerSylven/out/means.txt", check.names = FALSE)
pvals <- read.delim("/Users/ChristerSylven/out/pvalues.txt", check.names = FALSE)
#pval[pval==0] = 0.0009
#dim(means)
#dim(pvals)
# delete duplicates of L_R / R_L interactoin pairs
pvals0 <- pvals[!duplicated(pvals[,1]),]
#pvals1 <- pvals[unique(pvals[,1]),]
means0 <- means[!duplicated(means[,1]),]
#dim(pvals0)
#dim(means0)
#add L_R / R_L interaction pairs as rowname
rownames(pvals0) <- pvals0[,2]
rownames(means0) <- means0[,2]
#pvals0[1:2,1:10]
#means0[1:2,1:10]
# make dataframe w L_R R_L interaction pairs, secreted, mol 1 rec, mol 2 rec, integrin
pvals_LR <- pvals0[,c(7:9,11)]
pvals_LR <- cbind(rownames(pvals_LR), pvals_LR)
colnames(pvals_LR)[1] <- 'L_R0'
# to find number of interactions and number of integrin interactions for each muscle subtype and subset pvals0 below and get table on number of integrins
is_integrin <- function(cell_interactions) {
pvals2 <- pvals0[,cell_interactions]
# keep L_R R_L interaction pair where any value in row has p value < 0.01
pvals3 <- pvals2[apply(pvals2, 1, function(x) any(-log10(x) > 2)), ]
t2<- sum(pvals_LR[rownames(pvals3),]$is_integrin == 'False')
t3 <- sum(pvals_LR[rownames(pvals3),]$is_integrin == 'True')
t1 <- t2 + t3
# t[1,1] <- t1
t[1,1] <- t2
t[1,2] <- t3
t <<- t # makes dataframe global
}
plot_interactions_integrins <- function(Number) {
ordning <- c("Trabecular", "Compact", "Purkinje", "High G2M", "Exosome-related", 'OFT')
muscle_muscle <- data.frame(
Type = rep(c("nonintegrin", "integrin"), each = 6),
Cell_type = rep((ordning), 2),
Number = Number)
ggplot(muscle_muscle, aes(x = Cell_type, y = Number)) +
geom_col(aes(fill = Type), width = 0.7) +
scale_fill_manual(values = c("#0073C2FF", "#EFC000FF")) +
scale_x_discrete(limits = ordning, labels = ordning) +
ylim(0, 275) +
theme(
axis.text.x = element_text(angle = 90, hjust = 1, vjust = 0, size = 32),
axis.text.y = element_text(size = 32), legend.text = element_text(size = 32),
axis.title.x=element_blank(),
axis.title.y=element_blank(),
legend.title = element_text(size=24)
# axis.text.x=element_blank(),
# legend.position="none"
)
}
#theme(axis.line=element_blank(),
# axis.text.x=element_blank(),
#axis.text.y=element_blank(),
#axis.ticks=element_blank(),
#axis.title.x=element_blank(),
#axis.title.y=element_blank(),
Muscle_Muscle <- data.frame(
M0 = c( 'M1|M0', 'M2|M0', 'M4|M0', 'M6|M0', 'M7|M0', 'M0|M1', 'M0|M2', 'M0|M4', 'M0|M6', 'M0|M7'),
M1 = c( 'M0|M1', 'M2|M1', 'M4|M1', 'M6|M1', 'M7|M1', 'M1|M0', 'M1|M2', 'M1|M4', 'M1|M6', 'M1|M7'),
M2 = c( 'M0|M2', 'M2|M2', 'M4|M2', 'M6|M2', 'M7|M2', 'M2|M0', 'M2|M1', 'M2|M4', 'M2|M6', 'M2|M7'),
M4 = c( 'M0|M4', 'M1|M4', 'M2|M4', 'M6|M4', 'M7|M4', 'M4|M0', 'M4|M1', 'M4|M2', 'M4|M6', 'M4|M7'),
M6 = c( 'M0|M6', 'M1|M6', 'M2|M6', 'M4|M6', 'M7|M6', 'M6|M0', 'M6|M1', 'M6|M2', 'M6|M4', 'M6|M7'),
M7 =c( 'M0|M7', 'M1|M7', 'M2|M7', 'M4|M7', 'M6|M7', 'M7|M0', 'M7|M1', 'M7|M2', 'M7|M4', 'M7|M6')
)
t <- data.frame()
T <- data.frame()
for (i in 1:6) {
is_integrin( Muscle_Muscle[,i])
T <- rbind(T, t)
}
colnames(T) <- c('non_integrins', 'integrins')
T$clusters <- c('M0', 'M1', 'M2', 'M4', 'M6', 'M7')
#T
T2 <- melt(as.data.table(T), id = 'clusters')
#plot_interactions_integrins(T2$value)
ordning <- c("Trabecular", "Compact", "Purkinje", "High G2M", "Exosome-related", 'OFT')
muscle_muscle <- data.frame(
Type = rep(c("nonintegrin", "integrin"), each = 6),
Cell_type = rep((ordning), 2),
Number = T2$value)
p1 <- ggplot(muscle_muscle, aes(x = Cell_type, y = Number)) +
geom_col(aes(fill = Type), width = 0.7) +
scale_fill_manual(values = c("#0073C2FF", "#EFC000FF")) +
scale_x_discrete(limits = ordning, labels = ordning) +
ylim(0, 275) +
theme(
axis.text.x = element_blank(),
axis.text.y = element_text(size = 32), legend.text = element_text(size = 24),
axis.title.x=element_blank(),
axis.title.y=element_blank(),
legend.title = element_text(size=24),
# axis.text.x=element_blank(),
legend.position="none"
)
# Total /Integrin to non CM
# M0 -> M7
Muscle <- data.frame(
M0 = c('0|M0', '2|M0', '3|M0', '4|M0', '6|M0', '7|M0', '9|M0', '10|M0', '12|M0', '15|M0', '16|M0', 'M0|0', 'M0|2', 'M0|2', 'M0|4', 'M0|6', 'M0|7', 'M0|9', 'M0|10', 'M0|12', 'M0|15', 'M0|16'),
M1 = c('0|M1', '2|M1', '3|M1', '4|M1', '6|M1', '7|M1', '9|M1', '10|M1', '12|M1', '15|M1', '16|M1', 'M1|0', 'M1|2', 'M1|2', 'M1|4', 'M1|6', 'M1|7', 'M1|9', 'M1|10', 'M1|12', 'M1|15', 'M1|16'),
M2 = c('0|M2', '2|M2', '3|M2', '4|M2', '6|M2', '7|M2', '9|M2', '10|M2', '12|M2', '15|M2', '16|M2', 'M2|0', 'M2|2', 'M2|2', 'M2|4', 'M2|6', 'M2|7', 'M2|9', 'M2|10', 'M2|12', 'M2|15', 'M2|16'),
M4 = c('0|M4', '2|M4', '3|M4', '4|M4', '6|M4', '7|M4', '9|M4', '10|M4', '12|M4', '15|M4', '16|M4', 'M4|0', 'M4|2', 'M4|2', 'M4|4', 'M4|6', 'M4|7', 'M4|9', 'M4|10', 'M4|12', 'M4|15', 'M4|16'),
M6 = c('0|M6', '2|M6', '3|M6', '4|M6', '6|M6', '7|M6', '9|M6', '10|M6', '12|M6', '15|M6', '16|M6', 'M6|0', 'M6|2', 'M6|2', 'M6|4', 'M6|6', 'M6|7', 'M6|9', 'M6|10', 'M6|12', 'M6|15', 'M6|16'),
M7 =c('0|M7', '2|M7', '3|M7', '4|M7', '6|M7', '7|M7', '9|M7', '10|M7', '12|M7', '15|M7', '16|M7', 'M7|0', 'M7|2', 'M7|2', 'M7|4', 'M7|6', 'M7|7', 'M7|9', 'M7|10', 'M7|12', 'M7|15', 'M7|16')
)
t <- data.frame()
T <- data.frame()
for (i in 1:6) {
is_integrin( Muscle[,i])
T <- rbind(T, t)
}
colnames(T) <- c('non-integrins', 'integrins')
T$clusters <- c('M0', 'M1', 'M2', 'M4', 'M6', 'M7')
#T
T2 <- melt(as.data.table(T), id = 'clusters')
#plot_interactions_integrins(T2$value)
muscle_muscle <- data.frame(
Type = rep(c("nonintegrin", "integrin"), each = 6),
Cell_type = rep((ordning), 2),
Number = T2$value)
p2 <- ggplot(muscle_muscle, aes(x = Cell_type, y = Number)) +
geom_col(aes(fill = Type), width = 0.7) +
scale_fill_manual(values = c("#0073C2FF", "#EFC000FF")) +
scale_x_discrete(limits = ordning, labels = ordning) +
ylim(0, 275) +
theme(
#axis.text.x = element_text(angle = 90, hjust = 1, vjust = 0, size = 32),
axis.text.x = element_text(angle = 90, hjust = 1, vjust = 0, size = 24),
axis.text.y = element_text(size = 24), legend.text = element_text(size = 24),
axis.title.x=element_blank(),
axis.title.y=element_blank(),
legend.title = element_text(size=24),
# axis.text.x=element_blank(),
legend.position="bottom"
)
p1/p2To find Ligand_receptor dotplot according to cellphone DB and then edit results for
selected top LR pairs
rearrange so that pairs are shown as L-R and not R_L and make dotplot.
Plot as out/cellphoneDB plot and decide which pairs are top pairs
If out/cellphoneDB only show limited interactions as for M6, M7 analysis and ordering can be made directly manually.
system(‘cellphonedb plot dot_plot –rows in/OFT_pvals2_Nov_rows.txt –columns in/OFT_Nov.txt –output-name OFT_Nov.png’)
OFT_Nov.png saved in out folder with only 27 interactions, no selection is needed
#,2, 7, 6, 0, 3, 15
OFT <- c('2|M7', '7|M7', '6|M7', '0|M7', '3|M7', '15|M7', 'M7|2', 'M7|7', 'M7|6', 'M7|0', 'M7|3', 'M7|15')
write.table(OFT, file = '/Users/ChristerSylven/in/OFT_Nov.txt', sep = '\t', quote = F, row.names = F, col.names = F)
pvals2 <- pvals0[,OFT]
# to minimze number of interactions, filtered to those with p<0.01
pvals3 <- pvals2[apply(pvals2, 1, function(x) any(-log10(x) > 2)), ]
means3 <- means0[rownames(pvals3),OFT]
dim(pvals3) # 27 12## [1] 27 12
write.table(rownames(pvals3), '/Users/ChristerSylven/in/OFT_pvals2_Nov_rows.txt', sep = '\t', quote = F, row.names = F, col.names = F)
OFT <- c('2|M7', '7|M7', '6|M7', '0|M7', '3|M7', '15|M7', 'M7|2', 'M7|7', 'M7|6', 'M7|0', 'M7|3', 'M7|15')
pvals4_calc <- function(name, LRmrna, LRcells )
{
pvals2 <- pvals0[,LRcells] #pval0 duplicates eliminated above
pvals3 <- pvals2[apply(pvals2, 1, function(x) any(-log10(x) > 2)), ]
pvals3 <- pvals3[LRmrna,]
means3 <- means0[rownames(pvals3), LRcells]
pvals_LR_rows <- pvals_LR[rownames(pvals3),]
pvals4 <- cbind(pvals3,pvals_LR_rows)
means4 <- cbind(means3,pvals_LR_rows)
writexl::write_xlsx(pvals4, paste0("pvals4_d_", name, ".xlsx"))
writexl::write_xlsx(means4, paste0("means4_d_", name, ".xlsx"))
}
pvals4_calc('OFT', rownames(pvals3), OFT) # 27 pairs of which 12 are integrinsEdit in excel. Change RL to LR into L_R column. Change values accordingly.
For those changed correct mean molecule1/molecule2 in excel to molecule2/molecule1 (1/value)
Save as b for both corrected pvals and means files with all pairs as L_R.
pvals4b <- read_xlsx(paste0(path1, 'pvals4_d_OFT_b.xlsx'))
means4b <- read_xlsx(paste0(path1, 'means4_d_OFT_b.xlsx'))
# edit which pairs and order
#namn_OFT <- c(8:10,12,13:23,27, 24:25,28:30)
# function for manual dotplot, order and edited L_R pairs given by namn, default all LR pairs in 4b files
dotplot_LR <- function(pvals4b, means4b, L_Rs, length) # length <- length(OFT) etc
{
L_R <- length + 6
pvals4b <- as.data.frame(pvals4b, row.names=NULL)
means4b <- as.data.frame(means4b)
rownames(pvals4b) <- pvals4b[,L_R]
pvals4b <- pvals4b[,c(L_R, 1:length)]
pvals4b[pvals4b==0] = 0.0009
means4b <- means4b[,c(L_R, 1:length)]
means4b[means4b==0] = 1
pvals4b$L_R <- factor(pvals4b$L_R, levels = pvals4b$L_R )
means4b$L_R <- factor(means4b$L_R, levels = means4b$L_R )
pvals5 <- reshape2::melt(pvals4b, id=c('L_R'))
means5 <- reshape2::melt(means4b, id = c('L_R'))
level1 <- pvals5[1:dim(pvals4b)[1],1]
plot.data = cbind(pvals5,log2(means5[,3]))
colnames(plot.data) = c('pair', 'clusters', 'pvalue', 'mean')
#my_palette <- colorRampPalette(c("black", "blue", 'grey',"yellow", "red"), alpha=TRUE)(n=399)
my_palette <- colorRampPalette(c("black", "yellow","red"), alpha=TRUE)(n=399)
p <- ggplot2::ggplot(plot.data,aes(x=clusters,y=pair)) +
geom_point(aes(size=-log10(pvalue),color=mean)) +
scale_color_gradientn('Log2 mean (Ligand_receptor)', colors=my_palette) +
#scale_y_discrete(limits = level1, labels = level1)+
theme_bw() +
theme(panel.grid.minor = element_blank(),
panel.grid.major = element_blank(),
axis.text=element_text(size=12, colour = "black"),
axis.text.x = element_text(angle = 90, face = 'italic', hjust = 1, vjust=1),
axis.text.y = element_text(size=12, colour = "black"),
axis.title=element_blank(),
panel.border = element_rect(size = 0.7, linetype = "solid", colour = "black"),
legend.position="top",
plot.margin = unit(c(0.1,1,0,2), "cm")) +
scale_x_discrete(
labels = L_Rs) +
coord_flip()
#ggsave("/Users/ChristerSylven/Desktop/CMC_ms/GitHub_code/non_Cond_SAN_CM_L_R.pdf", p, height = 10, width = 12, device = "pdf")
p
}
L_Rs_OFT <- c("Epicardium_d~OFT",
"Fb_smooth_m~OFT",
"Fb_l_vascular~OFT",
"Endothelium~OFT",
"Fb_skeleton~OFT",
"Schwann~OFT",
"OFT~Epicardium_d",
"OFT~Fb_smooth_m",
"OFT~Fb_l_vascular",
"OFT~Endothelium",
"OFT~Fb_skeleton",
"OFT~Schwann")
dotplot_LR(pvals4b, means4b, L_Rs_OFT, length(OFT))non-integrin L-R interactions between Schwann progenitor cells and non CM celltypes colocalized in the outflow tract
Schwann 15 Epicardium_d 2, Fb_smooth-m 10, Fb_l_vascular 6, Endothelium 9, Fb_skeleton 3 All non-integrins interactions are shown.
Schwann <- c('2|15', '10|15', '6|15', '9|15', '3|15', '15|2', '15|10', '15|6', '15|9', '15|3')
pvals2 <- pvals0[,Schwann] #153
# only non-integrins
pvals0_1 <- pvals0[which( pvals0$is_integrin == 'False'),]
pvals2 <- pvals0_1[,Schwann]
# to minimze number of interactions, filtered to those with p<0.01
pvals3 <- pvals2[apply(pvals2, 1, function(x) any(-log10(x) > 2)), ]
means3 <- means0[rownames(pvals3),Schwann]
dim(pvals3) # 61 10## [1] 61 10
#add characteristics of each L-R
pvals_LR_Schwann <-pvals_LR[rownames(pvals3),]
pvals4 <- cbind(pvals3,pvals_LR_Schwann)
means4 <- cbind(means3,pvals_LR_Schwann)
#Schwann_non_integrin_rows2 L-R pairs ordered manually in Excel intocategories
Schwann_non_integrin_rows2 <- read.delim(paste0(path1,'data/Schwann_non_integrin_rows2.txt'), header = F)
#L-R pairs in Excel manually reordered into categories and saved as xlxx file
pvals4_reorder <- pvals4[match(Schwann_non_integrin_rows2[,1], rownames(pvals4)),]
means4_reorder <- means4[match(Schwann_non_integrin_rows2[,1], rownames(means4)),]
writexl::write_xlsx(pvals4_reorder, paste0(path1,"/pvals4_", 'Schwann', ".xlsx"))
writexl::write_xlsx(means4_reorder, paste0(path1,"/means4_", 'Schwann', ".xlsx"))
# edit in excel R_L to L_R, add column L_R with correct L_R names, invert transformed values in means, save as:pvals4_Schwann_d.xlsx, means4_Schwann_d.xlsx
L_Rs_Schwann <- c(
"Epicardium_d~Schwann",
"Fb_smooth_m~Schwann",
"Fb_l_vascular~Schwann",
"Endothelium~Schwann",
"Fb_skeleton~Schwann",
"Schwann~Epicardium_d",
"Schwann~Fb_smooth_m",
"Schwann~Fb_l_vascular",
"Schwann~Endothelium",
"Schwann~Fb_skeleton"
)
pvals4b <- read_xlsx(paste0(path1,'pvals4_Schwann_d.xlsx'))
means4b <- read_xlsx(paste0(path1,'means4_Schwann_d.xlsx'))
dotplot_LR(pvals4b, means4b, L_Rs_Schwann, length(Schwann))EXO <- c('2|M6', '7|M6', '6|M6', '0|M6', '3|M6', '15|M6', 'M6|2', 'M6|7', 'M6|6', 'M6|0', 'M6|3', 'M6|15')
write.table(EXO, file = './in/EXO_Nov.txt', sep = '\t', quote = F, row.names = F, col.names = F)
pvals2 <- pvals0[,EXO]
# to minimze number of interactions, filtered to those with p<0.01
pvals3 <- pvals2[apply(pvals2, 1, function(x) any(-log10(x) > 2)), ]
means3 <- means0[rownames(pvals3),EXO]
dim(pvals3) # 31 12## [1] 32 12
write.table(rownames(pvals3), './in/EXO_pvals2_Nov_rows.txt', sep = '\t', quote = F, row.names = F, col.names = F)
pvals4_calc('EXO', rownames(pvals3), EXO) # 32 pairs of which 11 are integrinssystem(‘cellphonedb plot dot_plot –rows in/EXO_pvals2_Nov_rows.txt –columns in/EXO_Nov.txt –output-name EXO_Nov.png’)
# edit in excel. Change RL to LR into L_R column. Change values accordingly. Save as b for both files.
# Correct mean molecule1/molecule2 in excel with to molecule2/molecule1 (1/value)
#
#corrected files with all pairs as L_R saved as b
pvals4b <- read_xlsx(paste0(path1, 'pvals4_d_EXO_b.xlsx'))
means4b <- read_xlsx(paste0(path1, 'means4_d_EXO_b.xlsx'))
L_Rs_EXO <- c("Epicardium_d~EXO",
"Fb_smooth_m~EXO",
"Fb_l_vascular~EXO",
"Endothelium~EXO",
"Fb_skeleton~EXO",
"Schwann~EXO",
"EXO~Epicardium_d",
"EXO~Fb_smooth_m",
"EXO~Fb_l_vascular",
"EXO~Endothelium",
"EXO~Fb_skeleton",
"EXO~Schwann")
# Edited L_R order:
# namn_exo <- c(2:19,23:25,22, 28:29,20:21,27,30,32:33,26,31)
dotplot_LR(pvals4b, means4b, L_Rs_EXO, length(EXO))# calculate number of L_Rs and integrin L_Rs in Conduit CMs
#
Conduit_lig<- c('M5|0', 'M5|2', 'M5|6', 'M5|7', 'M5|12', 'M5|15')
pvals2 <- pvals0[,Conduit_lig]
pvals3 <- pvals2[apply(pvals2, 1, function(x) any(-log10(x) > 2)), ]
#dim(pvals3) # 102
LR_cond <- dim(pvals3)[1]
pvals0_1 <- pvals0[which( pvals0$is_integrin == 'False'),]
pvals2_1 <- pvals0_1[,Conduit_lig]
pvals3_1 <- pvals2_1[apply(pvals2_1, 1, function(x) any(-log10(x) > 2)), ]
#dim(pvals3_1) # 70
LR_cond_non_integrin <- dim(pvals3_1)[1]
Conduit_rec<- c('0|M5', '2|M5', '6|M5', '7|M5', '12|M5', '15|M5')
pvals2 <- pvals0[,Conduit_rec]
pvals3 <- pvals2[apply(pvals2, 1, function(x) any(-log10(x) > 2)), ]
#dim(pvals3) # 124
LR_cond <- LR_cond + dim(pvals3)[1]
pvals0_1 <- pvals0[which( pvals0$is_integrin == 'False'),]
pvals2_1 <- pvals0_1[,Conduit_rec]
pvals3_1 <- pvals2_1[apply(pvals2_1, 1, function(x) any(-log10(x) > 2)), ]
#dim(pvals3_1)
LR_cond_non_integrin <- LR_cond_non_integrin + dim(pvals3_1)[1]
# VCAM1 Ig adhesion molecule viewed to have non integrin function and added to non integrins
cat('Number of Conduit CM total L-R interactions: ', LR_cond, ' and integrin interactions: ', (LR_cond - LR_cond_non_integrin - 2),
'\nVCAM1 Ig adhesion molecule viewed to have non integrin function and added to non integrins.',
'\nFigure shows L-R interactions for non integrins and for non collagen integrins.')## Number of Conduit CM total L-R interactions: 226 and integrin interactions: 72
## VCAM1 Ig adhesion molecule viewed to have non integrin function and added to non integrins.
## Figure shows L-R interactions for non integrins and for non collagen integrins.
# edited order of Conduit L-R and R-l receptor pairs. Non collagen integrins added (n=10)
COND_lig_rec <- read_xlsx(paste0(path1, 'data/COND_lig_rec.xlsx')) # Edited order of non collagen L_R pairs for Counduit (L-R) and Conduit (R_L)
CONDUIT_lig <- as.data.frame(na.omit(COND_lig_rec[,1]))
pvals4_calc('Cond_lig2', CONDUIT_lig[,1], Conduit_lig)
CONDUIT_rec <- as.data.frame(na.omit(COND_lig_rec[,2]))
pvals4_calc('Cond_rec2', CONDUIT_rec[,1], Conduit_rec)
# edit in excel. Change RL to LR into L_R column. Change values accordingly. Save as b for both pvalue4 and mean4 files.
# Correct mean molecule1/molecule2 in excel with to molecule2/molecule1 (1/value)
#
#corrected files with all pairs as L_R saved as b
#
L_Rs_nonCM_COND_lig <- c(
"Conduit~Endothelium",
"Conduit~Epicardium_d",
"Conduit~Fb_l_vascular",
"Conduit~Fb_smooth_m",
"Conduit~Epicardium",
"Conduit-Schwann")
L_Rs_nonCM_COND_rec <- c(
"Endothelium~Conduit",
"Epicardium_d~Conduit",
"Fb_l_vascular~Conduit",
"Fb_smooth_m~Conduit",
"Epicardium~Conduit",
"Schwann~Conduit")pvals4b <- read_xlsx(paste0(path1, 'pvals4_d_COND_lig_b.xlsx'))#87
means4b <- read_xlsx(paste0(path1, 'means4_d_COND_lig_b.xlsx'))
d1 <- dotplot_LR(pvals4b, means4b,L_Rs_nonCM_COND_lig, length(Conduit_lig))
pvals4b <- read_xlsx(paste0(path1, 'pvals4_d_COND_rec_b.xlsx')) #75 87+75 = 162 CONDUIT_lig (66) + CONDUIT_rec (80) = 146 + 16? 11 non coll integrins
means4b <- read_xlsx(paste0(path1, 'means4_d_COND_rec_b.xlsx'))
d2 <- dotplot_LR(pvals4b, means4b,L_Rs_nonCM_COND_rec, length(Conduit_rec))
d1 / d2sessionInfo() ## R version 4.2.0 (2022-04-22)
## Platform: x86_64-apple-darwin17.0 (64-bit)
## Running under: macOS Big Sur/Monterey 10.16
##
## Matrix products: default
## BLAS: /Library/Frameworks/R.framework/Versions/4.2/Resources/lib/libRblas.0.dylib
## LAPACK: /Library/Frameworks/R.framework/Versions/4.2/Resources/lib/libRlapack.dylib
##
## locale:
## [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
##
## attached base packages:
## [1] stats graphics grDevices utils datasets methods base
##
## other attached packages:
## [1] writexl_1.4.0 xlsx_0.6.5 readxl_1.4.0 openxlsx_4.2.5
## [5] reshape2_1.4.4 data.table_1.14.2 dplyr_1.0.9 sp_1.5-0
## [9] SeuratObject_4.1.0 Seurat_4.1.1 biomaRt_2.52.0 ggplot2_3.3.6
## [13] igraph_1.3.4
##
## loaded via a namespace (and not attached):
## [1] BiocFileCache_2.4.0 plyr_1.8.7 lazyeval_0.2.2
## [4] splines_4.2.0 listenv_0.8.0 scattermore_0.8
## [7] GenomeInfoDb_1.32.2 digest_0.6.29 htmltools_0.5.3
## [10] fansi_1.0.3 magrittr_2.0.3 memoise_2.0.1
## [13] tensor_1.5 cluster_2.1.3 ROCR_1.0-11
## [16] globals_0.15.1 Biostrings_2.64.0 matrixStats_0.62.0
## [19] spatstat.sparse_2.1-1 prettyunits_1.1.1 colorspace_2.0-3
## [22] blob_1.2.3 rappdirs_0.3.3 ggrepel_0.9.1
## [25] xfun_0.31 crayon_1.5.1 RCurl_1.98-1.7
## [28] jsonlite_1.8.0 progressr_0.10.1 spatstat.data_2.2-0
## [31] survival_3.3-1 zoo_1.8-10 glue_1.6.2
## [34] polyclip_1.10-0 gtable_0.3.0 zlibbioc_1.42.0
## [37] XVector_0.36.0 leiden_0.4.2 future.apply_1.9.0
## [40] BiocGenerics_0.42.0 abind_1.4-5 scales_1.2.0
## [43] DBI_1.1.3 spatstat.random_2.2-0 miniUI_0.1.1.1
## [46] Rcpp_1.0.9 viridisLite_0.4.0 xtable_1.8-4
## [49] progress_1.2.2 reticulate_1.25-9000 spatstat.core_2.4-4
## [52] bit_4.0.4 stats4_4.2.0 htmlwidgets_1.5.4
## [55] httr_1.4.3 RColorBrewer_1.1-3 ellipsis_0.3.2
## [58] ica_1.0-3 farver_2.1.1 rJava_1.0-6
## [61] pkgconfig_2.0.3 XML_3.99-0.10 uwot_0.1.11
## [64] deldir_1.0-6 sass_0.4.2 dbplyr_2.2.1
## [67] utf8_1.2.2 labeling_0.4.2 tidyselect_1.1.2
## [70] rlang_1.0.4 later_1.3.0 AnnotationDbi_1.58.0
## [73] cellranger_1.1.0 munsell_0.5.0 tools_4.2.0
## [76] cachem_1.0.6 cli_3.3.0 generics_0.1.3
## [79] RSQLite_2.2.15 ggridges_0.5.3 evaluate_0.15
## [82] stringr_1.4.0 fastmap_1.1.0 goftest_1.2-3
## [85] yaml_2.3.5 knitr_1.39 bit64_4.0.5
## [88] fitdistrplus_1.1-8 zip_2.2.0 purrr_0.3.4
## [91] RANN_2.6.1 KEGGREST_1.36.3 nlme_3.1-158
## [94] pbapply_1.5-0 future_1.27.0 mime_0.12
## [97] xml2_1.3.3 compiler_4.2.0 rstudioapi_0.13
## [100] plotly_4.10.0 filelock_1.0.2 curl_4.3.2
## [103] png_0.1-7 spatstat.utils_2.3-1 tibble_3.1.8
## [106] bslib_0.4.0 stringi_1.7.8 highr_0.9
## [109] rgeos_0.5-9 lattice_0.20-45 Matrix_1.4-1
## [112] vctrs_0.4.1 pillar_1.8.0 lifecycle_1.0.1
## [115] spatstat.geom_2.4-0 lmtest_0.9-40 jquerylib_0.1.4
## [118] RcppAnnoy_0.0.19 cowplot_1.1.1 bitops_1.0-7
## [121] irlba_2.3.5 httpuv_1.6.5 patchwork_1.1.1
## [124] R6_2.5.1 promises_1.2.0.1 KernSmooth_2.23-20
## [127] gridExtra_2.3 IRanges_2.30.0 parallelly_1.32.1
## [130] codetools_0.2-18 MASS_7.3-58 assertthat_0.2.1
## [133] xlsxjars_0.6.1 withr_2.5.0 sctransform_0.3.3
## [136] S4Vectors_0.34.0 GenomeInfoDbData_1.2.8 mgcv_1.8-40
## [139] parallel_4.2.0 hms_1.1.1 rpart_4.1.16
## [142] grid_4.2.0 tidyr_1.2.0 rmarkdown_2.14
## [145] Rtsne_0.16 Biobase_2.56.0 shiny_1.7.2